1 Background

Considering the prevalence of obesity (1), its associated fracture risk [(2)] and that these patients are especially prone to develop loading associated musculoskeletal injuries [(5)], there is a need to develop mechanical loading prediction equations that are accurate for overweight and obese subjects in order to precisely determine and monitor exercise-associated mechanical loading in these patients. Therefore, the purpose of this study was to develop prediction equations based on accelerometer data able to accurately predict ground reaction forces (GRF) and loading rate (LR) on a broad range of body masses, from normal weight to severely obese subjects, setting thereby, a base for objective prescription and monitoring of exercise mechanical loads.

2 Load data

library(knitr)
library(broman)
library(tidyverse)
library(here)
library(nlme)
library(piecewiseSEM)
library(ez)
source(here("R", "BMI_categories.R"))
source(here("R", "cross_validate_mixed_model.R"))
source(here("R", "accuracy_indices.R"))
source(here("R", "get_BA_plot.R"))

Initially, 4 files needed to be loaded, containing data from each of the accelerometer placements utilized in the study (lower back and hip), and for each of the mechanical loading variables (GRF and LR).

back <- read_csv(here("data", "back_sec.csv")) %>% 
  select(ID, speed, body_mass, height, BMI, sex, age, pVGRF_N, pVACC_g, pRGRF_N, pRACC_g) %>% 
  filter(speed <= 6)

hip <- read_csv(here("data", "hip_sec.csv")) %>% 
  select(ID, speed, body_mass, height, BMI, sex, age, pVGRF_N, pVACC_g, pRGRF_N, pRACC_g) %>% 
  filter(speed <= 6)

anthopometric <- read_csv(here("data", "back_sec.csv")) %>% 
  select(ID, speed, height, BMI, sex, age)

back_LR <- read_csv(here("data", "max_rates_back_IMU.csv")) %>% 
  left_join(anthopometric, by = c("ID", "speed")) %>% 
  BMI_categories() %>% 
  select(
    ID, speed, body_mass, height, BMI, BMI_cat, sex, age, pVLR_Ns, pVATR_gs, pRLR_Ns, pRATR_gs
  )

hip_LR <- read_csv(here("data", "max_rates_hip_IMU.csv")) %>% 
  left_join(anthopometric, by = c("ID", "speed")) %>% 
  BMI_categories() %>% 
  select(
    ID, speed, body_mass, height, BMI, BMI_cat, sex, age, pVLR_Ns, pVATR_gs, pRLR_Ns, pRATR_gs
  )

A piece of a data frame used in all subsequent analysis can be seen below.

# For GRF
back
## # A tibble: 300 x 11
##       ID speed body_mass height   BMI sex     age pVGRF_N pVACC_g pRGRF_N
##    <dbl> <dbl>     <dbl>  <dbl> <dbl> <chr> <dbl>   <dbl>   <dbl>   <dbl>
##  1     3     2      85.6    171  29.3 M        50   1023.    1.26   1030.
##  2     3     3      85.6    171  29.3 M        50   1118.    1.51   1133.
##  3     3     4      85.6    171  29.3 M        50   1202.    1.61   1230.
##  4     3     5      85.6    171  29.3 M        50   1200.    1.63   1246.
##  5     3     6      85.6    171  29.3 M        50   1232.    1.73   1294.
##  6     4     2      61.1    162  23.3 F        45    671.    1.16    678.
##  7     4     3      61.1    162  23.3 F        45    702.    1.23    716.
##  8     4     4      61.1    162  23.3 F        45    756.    1.35    770.
##  9     4     5      61.1    162  23.3 F        45    809.    1.50    828.
## 10     4     6      61.1    162  23.3 F        45    912.    1.67    926.
## # … with 290 more rows, and 1 more variable: pRACC_g <dbl>
# For LR
back_LR
## # A tibble: 271 x 12
##       ID speed body_mass height   BMI BMI_cat sex     age pVLR_Ns pVATR_gs
##    <dbl> <dbl>     <dbl>  <dbl> <dbl> <fct>   <chr> <dbl>   <dbl>    <dbl>
##  1     3     2      86.5    171  29.3 overwe… M        50   4367.     7.95
##  2     3     3      86.5    171  29.3 overwe… M        50   7895.    13.3 
##  3     3     4      86.5    171  29.3 overwe… M        50  10585.    14.1 
##  4     3     5      86.5    171  29.3 overwe… M        50  12798.    15.1 
##  5     3     6      86.5    171  29.3 overwe… M        50  13556.    20.2 
##  6     4     2      61.1    162  23.3 normal… F        45   3650.     4.12
##  7     4     3      61.1    162  23.3 normal… F        45   5049.     5.42
##  8     4     4      61.1    162  23.3 normal… F        45   6872.    11.3 
##  9     4     5      61.1    162  23.3 normal… F        45   7992.    17.4 
## 10     4     6      61.1    162  23.3 normal… F        45  10800.    15.3 
## # … with 261 more rows, and 2 more variables: pRLR_Ns <dbl>,
## #   pRATR_gs <dbl>

3 Mechanical loading X acceleration relationships

We first determined the relationship between peak acceleration (pACC) and peak GRF (pGRF), and between peak ACC acceleration transient rate (pATR) and peak LR (pLR) obtained during the incremental walking speeds used in our experimental protocol. This was performed in order to verify the consistency of these relationships for increasing pACC and pATR values and the ability of the accelerometer to discriminate differences in pGRF and pLR between subjects in different body mass indices (BMI) classes. A scatterplot with these relationships for all two accelerometer placements tested and for both the resultant and its vertical component is depicted in Figure 1. Generally, as expected, pACC recorded by accelerometers in all placements, showed a linear increase as a function of pGRF increases. Also, the recorded ACC were shown to be able to provide a good discrimination between different BMI classes, as for the same registered pACC the pGRF tended to be consistently higher for subjects in higher BMI classes. As was also expected, the pATR recorded by the accelerometers demonstrated a linear increase as the pLR increases. Contrarily to what was observed on the relationship between pACC and pGRF, the pATR x pLR plot does not provide a good discrimination between different BMI classes. Moreover, the values at low magnitudes showed a high degree of clustering, whereas the values at higher magnitudes tended to a greater dispersion.

source(here("figs", "fig1.R"))

GRF_ACC_plot_grid

Code to generate these plots can be seen here.

4 Ground reaction force prediction

4.1 Newton’s second law of motion

In a first attempt to predict the pGRF using accelerometry data, Newton’s second law of motion (force = mass · acceleration) was applied to predict peak resultant GRF (pRGRF) from peak resultant ACC (pRACC) derived from each accelerometer placement. As both accelerometers were located at waistline, both were close to the body center of mass, and therefore, their recorded ACC values could represent the ACC of the body center of mass.

The code usede to make these predictions can be seen below.

# Back
back_2nd_law <- back %>% 
  mutate(
    pRACC_ms2 = pRACC_g * 9.81,
    pRGRF_N_predicted = body_mass * pRACC_ms2
  )

# Hip
hip_2nd_law <- hip %>% 
  mutate(
    pRACC_ms2 = pRACC_g * 9.81,
    pRGRF_N_predicted = body_mass * pRACC_ms2
  )

To evaluate these predictions, some accuracy indices were computed. These indices were the mean absolute error (MAE), mean absolute percent error (MAPE) and root mean square error (RMSE), and they were computed with the accuracy_indices() function.

# Back
back_2nd_law_accuracy <- accuracy_indices(back_2nd_law, "pRGRF_N", "pRGRF_N_predicted")
back_2nd_law_accuracy
##   bias_mean bias_sd LoA_inf LoA_sup    MAE MAE_sd   RMSE bias_mean_p
## 1    -90.37  137.71 -360.28  179.54 109.86 122.68 164.52       -7.21
##   bias_sd_p LoA_inf_p LoA_sup_p MAPE MAPE_sd CV_RMSE
## 1      10.7    -28.18     13.76 9.21    9.03   14.06
# Hip
hip_2nd_law_accuracy <- accuracy_indices(hip_2nd_law, "pRGRF_N", "pRGRF_N_predicted")
hip_2nd_law_accuracy
##   bias_mean bias_sd LoA_inf LoA_sup   MAE MAE_sd   RMSE bias_mean_p
## 1   -224.65   191.6 -600.19  150.89 229.4 185.86 295.05      -18.31
##   bias_sd_p LoA_inf_p LoA_sup_p  MAPE MAPE_sd CV_RMSE
## 1     13.86    -45.48      8.86 18.85   13.11   25.25

4.2 Linear mixed models

Linear mixed models (LMM) were developed to predict peak resultant GRF (pRGRF) and peak vertical GRF (pVGRF). Distinct LMMs were developed with data from lowee back and hip accelerometer placement (back and hip data frames) using the lme() function of the nlme package. Covariance structure used was an autoregressive process of order 1 (correlation = corAR1) and maximum likelihood method was used for estimating parameters (method = "ML"). Predictors tested on pRGRF models were body mass and peak resultant acceleration (pRACC), while on pVGRF were body mass and vertical acceleration (pVACC). All of them were tested as fixed effects and have shown to be significant predictors (e.g.: fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass). Both random intercept and slopes were tested, but only the random intercept inclusion has showed models improvement (random = ~ 1 | ID). Linear, quadratic and cubic polynomial simulations were also tested, whereas the last one did not contribute significantly to the models. Final models were chosen according to -2 log-likelihood statistics. Traditional coefficient of determination (R2) was represented by conditional R2 (6) computed with rsquared() function of the piecewiseSEM package, that estimates the variance explained by the whole model (7).

Code to build the LMMs, as well as their summary and R2, can be found below.

# For resultant peak ground reaction force
back_res_LMM <- lme(
  fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g : body_mass,
  random = ~ 1 | ID,
  method = "ML",
  correlation = corAR1(),
  data = back
)
r2_back_res_LMM <- rsquared(back_res_LMM)

summary(back_res_LMM)
## Linear mixed-effects model fit by maximum likelihood
##  Data: back 
##        AIC      BIC    logLik
##   3185.997 3215.628 -1584.999
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.0184402 79.53301
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.8530994 
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)       -698.7031 123.44832 233 -5.659884   0e+00
## pRACC_g           1047.5129 122.40109 233  8.558036   0e+00
## I(pRACC_g^2)      -345.2605  36.10073 233 -9.563809   0e+00
## body_mass            3.8294   1.06086  62  3.609729   6e-04
## pRACC_g:body_mass    6.0219   0.61893 233  9.729449   0e+00
##  Correlation: 
##                   (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g           -0.896                      
## I(pRACC_g^2)       0.648 -0.892               
## body_mass         -0.636  0.268   0.143       
## pRACC_g:body_mass  0.565 -0.281  -0.167 -0.922
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.0955967 -0.6001009 -0.0930397  0.5792458  3.1083008 
## 
## Number of Observations: 300
## Number of Groups: 64
r2_back_res_LMM
##   Response   family     link method  Marginal Conditional
## 1  pRGRF_N gaussian identity   none 0.9361655   0.9361655
hip_res_LMM <- lme(
  fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g : body_mass,
  random = ~ 1 | ID,
  method = "ML",
  correlation = corAR1(),
  data = hip
)
r2_hip_res_LMM <- rsquared(hip_res_LMM)


summary(hip_res_LMM)
## Linear mixed-effects model fit by maximum likelihood
##  Data: hip 
##        AIC      BIC    logLik
##   3091.034 3120.637 -1537.517
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:  0.02080893  71.4172
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.8662087 
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass 
##                       Value Std.Error  DF    t-value p-value
## (Intercept)       -300.9909  89.71869 232  -3.354829   9e-04
## pRACC_g            522.6850  71.61994 232   7.298037   0e+00
## I(pRACC_g^2)      -171.5606  16.91100 232 -10.144910   0e+00
## body_mass            3.9596   0.82234  62   4.814986   0e+00
## pRACC_g:body_mass    5.3671   0.41930 232  12.799992   0e+00
##  Correlation: 
##                   (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g           -0.878                      
## I(pRACC_g^2)       0.588 -0.854               
## body_mass         -0.766  0.422   0.030       
## pRACC_g:body_mass  0.678 -0.471  -0.038 -0.891
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.54273847 -0.59285086  0.05857753  0.70507501  2.51256437 
## 
## Number of Observations: 299
## Number of Groups: 64
r2_hip_res_LMM
##   Response   family     link method  Marginal Conditional
## 1  pRGRF_N gaussian identity   none 0.9488446   0.9488446
# For vertical peak ground reaction force
back_vert_LMM <- lme(
  fixed = pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g : body_mass,
  random = ~ 1 | ID,
  method = "ML",
  correlation = corAR1(),
  data = back
)
r2_back_vert_LMM <- rsquared(back_vert_LMM)

summary(back_vert_LMM)
## Linear mixed-effects model fit by maximum likelihood
##  Data: back 
##        AIC      BIC    logLik
##   3234.998 3264.628 -1609.499
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:  0.02508531 102.5201
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.9079104 
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)       -776.8934 145.57197 233 -5.336834       0
## pVACC_g           1042.9052 146.96000 233  7.096524       0
## I(pVACC_g^2)      -336.2115  44.82805 233 -7.500024       0
## body_mass            6.2132   1.22304  62  5.080120       0
## pVACC_g:body_mass    5.0805   0.72861 233  6.972936       0
##  Correlation: 
##                   (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g           -0.889                      
## I(pVACC_g^2)       0.652 -0.894               
## body_mass         -0.688  0.327   0.066       
## pVACC_g:body_mass  0.587 -0.341  -0.101 -0.891
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2287736 -0.6177471 -0.1216603  0.4585545  3.7585643 
## 
## Number of Observations: 300
## Number of Groups: 64
r2_back_vert_LMM
##   Response   family     link method  Marginal Conditional
## 1  pVGRF_N gaussian identity   none 0.8935621   0.8935621
hip_vert_LMM <- lme(
  fixed = pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g : body_mass,
  random = ~ 1 | ID,
  method = "ML",
  correlation = corAR1(),
  data = hip
)
r2_hip_vert_LMM <- rsquared(hip_vert_LMM)

summary(hip_vert_LMM)
## Linear mixed-effects model fit by maximum likelihood
##  Data: hip 
##        AIC      BIC    logLik
##   3153.414 3183.017 -1568.707
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:  0.02202366 84.32584
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.8869442 
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)       -435.7365 122.69794 232 -3.551294   5e-04
## pVACC_g            586.6627 112.66171 232  5.207295   0e+00
## I(pVACC_g^2)      -188.9689  28.55314 232 -6.618146   0e+00
## body_mass            5.8047   0.97485  62  5.954511   0e+00
## pVACC_g:body_mass    4.9544   0.54067 232  9.163472   0e+00
##  Correlation: 
##                   (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g           -0.910                      
## I(pVACC_g^2)       0.715 -0.907               
## body_mass         -0.759  0.471  -0.125       
## pVACC_g:body_mass  0.688 -0.524   0.135 -0.888
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.95040403 -0.53289689 -0.04360304  0.45305786  4.09646433 
## 
## Number of Observations: 299
## Number of Groups: 64
r2_hip_vert_LMM
##   Response   family     link method  Marginal Conditional
## 1  pVGRF_N gaussian identity   none 0.9263472   0.9263473

In all models, conditional R2 values ranged from 0.89 to 0.95.

4.3 Validation analysis

Model validation was assessed by the leave-one-out cross-validation (LOOCV) method. This approach is recommended when a different sample is not available for cross-validation and it provides an insight on the model potential to predict outcomes in a new independent sample (8). For LOOCV each participant’s data was separated into a testing dataset (one at a time) with the remaining data being in the training dataset. New LMM, with the same outcomes and predictors as determined for the entire sample, were developed using the training dataset and then used to predict pGRF for the only participant in the testing dataset. This process was repeated for all participants (64 times). Data from testing dataset were used in the remaining statistical analysis.

LOOCV of the LMM was done with the cross_validate_mixed_model() function and is shown below.

# For resultant peak ground reaction force
fix_eff    <- pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass
rand_eff   <- ~ 1 | ID
# Back
LOOCV_back_res_LMM <- do.call(rbind, (lapply(unique(back$ID), cross_validate_mixed_model, df = back)))
# Hip
LOOCV_hip_res_LMM <- do.call(rbind, (lapply(unique(hip$ID), cross_validate_mixed_model, df = hip)))

# For vertical peak ground reaction force
fix_eff    <- pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass
rand_eff   <- ~ 1 | ID
# Back
LOOCV_back_vert_LMM <- do.call(rbind, (lapply(unique(back$ID), cross_validate_mixed_model, df = back)))
# Hip
LOOCV_hip_vert_LMM <- do.call(rbind, (lapply(unique(hip$ID), cross_validate_mixed_model, df = hip)))

4.3.1 Bland-Altman plots

Bland-Altman plots were used to examine the agreement between pGRF measured with force plates and those predicted through the regression equations. The difference of the actual and predicted pGRF was plotted against their mean. Bias was expressed as the mean of these differences and the limits of agreement were obtained using ±1.96 standard deviation of the mean between actual and predicted pGRF (9).

Bland-Altman plots can be seen below for all accelerometer placements. Panel A shows plots of pRGRF and panel B shows plots of pVGRF.

source(here("figs", "fig2.R"))

BA_plot_grid_1 <- plot_grid(
  back_res_BA_plot     + theme(legend.position = "none"),
  back_vert_BA_plot    + theme(legend.position = "none"),
  hip_res_BA_plot      + theme(legend.position = "none"),
  hip_vert_BA_plot     + theme(legend.position = "none"),
  labels = c("A", "B", "", ""),
  align  = "h", vjust = 1, label_size = 16,
  ncol   = 2, nrow = 2
)

legend <- get_legend(back_res_BA_plot)

BA_plot_grid <- plot_grid(BA_plot_grid_1, legend, ncol = 1, rel_heights = c(1, 0.1))

BA_plot_grid

Code to generate these plots can be seen here.


One-sample t-tests were run to check whether bias from each accelerometer placement, in both pRGRF and pVGRF, was significantly different than 0. These tests were run using the t.test() function of the base R with the argument mu = 0. No significant differences were found.

# For resultant peak ground reaction force
# Back
LOOCV_back_res_LMM$diff <- LOOCV_back_res_LMM$pRGRF_N - LOOCV_back_res_LMM$pRGRF_N_predicted
LOOCV_back_res_LMM$mean <- (LOOCV_back_res_LMM$pRGRF_N + LOOCV_back_res_LMM$pRGRF_N_predicted) / 2
t.test(LOOCV_back_res_LMM$diff, mu = 0)
## 
##  One Sample t-test
## 
## data:  LOOCV_back_res_LMM$diff
## t = 0.44993, df = 299, p-value = 0.6531
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -7.724802 12.303988
## sample estimates:
## mean of x 
##  2.289593
# Hip
LOOCV_hip_res_LMM$diff <- LOOCV_hip_res_LMM$pRGRF_N - LOOCV_hip_res_LMM$pRGRF_N_predicted
LOOCV_hip_res_LMM$mean <- (LOOCV_hip_res_LMM$pRGRF_N + LOOCV_hip_res_LMM$pRGRF_N_predicted) / 2
t.test(LOOCV_hip_res_LMM$diff, mu = 0)
## 
##  One Sample t-test
## 
## data:  LOOCV_hip_res_LMM$diff
## t = 0.36236, df = 298, p-value = 0.7173
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -7.774785 11.284133
## sample estimates:
## mean of x 
##  1.754674
# For vertical peak ground reaction force
# Back
LOOCV_back_vert_LMM$diff <- LOOCV_back_vert_LMM$pVGRF_N - LOOCV_back_vert_LMM$pVGRF_N_predicted
LOOCV_back_vert_LMM$mean <- (LOOCV_back_vert_LMM$pVGRF_N + LOOCV_back_vert_LMM$pVGRF_N_predicted) / 2
t.test(LOOCV_back_vert_LMM$diff, mu = 0)
## 
##  One Sample t-test
## 
## data:  LOOCV_back_vert_LMM$diff
## t = -0.10563, df = 299, p-value = 0.9159
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -12.19819  10.95537
## sample estimates:
##  mean of x 
## -0.6214066
# Hip
LOOCV_hip_vert_LMM$diff <- LOOCV_hip_vert_LMM$pVGRF_N - LOOCV_hip_vert_LMM$pVGRF_N_predicted
LOOCV_hip_vert_LMM$mean <- (LOOCV_hip_vert_LMM$pVGRF_N + LOOCV_hip_vert_LMM$pVGRF_N_predicted) / 2
t.test(LOOCV_hip_vert_LMM$diff, mu = 0)
## 
##  One Sample t-test
## 
## data:  LOOCV_hip_vert_LMM$diff
## t = 0.0087888, df = 298, p-value = 0.993
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -9.956701 10.046032
## sample estimates:
##  mean of x 
## 0.04466578

Also, linear regressions were applied to identify if there was any proportional bias, that is, if bias was related with the magnitude of the mean between measured and predicted pGRF (10). Linear regressions were run using the lm() function of the base R.

Results showed a significant proportional bias (p < 0.05), however, with a low magnitude (highest R2 = 0.032). These results suggest that despite there is a trend for underestimation at increasingly higher pGRF values, the magnitude of this effect is neglectable.

# For resultant peak ground reaction force
# Back
back_res_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_res_LMM)
summary(back_res_BA_plot_LR)
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_back_res_LMM)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -267.995  -55.071   -2.505   52.450  294.243 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -47.77144   19.69289  -2.426  0.01587 * 
## mean          0.04282    0.01628   2.630  0.00899 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 87.28 on 298 degrees of freedom
## Multiple R-squared:  0.02268,    Adjusted R-squared:  0.0194 
## F-statistic: 6.915 on 1 and 298 DF,  p-value: 0.008991
# Hip
hip_res_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_hip_res_LMM)
summary(hip_res_BA_plot_LR)
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_hip_res_LMM)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -268.721  -53.909    5.937   51.546  235.762 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -40.38734   18.75117  -2.154   0.0321 *
## mean          0.03610    0.01552   2.325   0.0207 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 83.12 on 297 degrees of freedom
## Multiple R-squared:  0.01788,    Adjusted R-squared:  0.01457 
## F-statistic: 5.406 on 1 and 297 DF,  p-value: 0.02074
# For vertical peak ground reaction force
# Back
back_vert_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_vert_LMM)
summary(back_vert_BA_plot_LR)
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_back_vert_LMM)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -342.85  -64.14   -4.92   51.55  354.11 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -46.47386   23.16029  -2.007   0.0457 *
## mean          0.03998    0.01954   2.046   0.0416 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 101.4 on 298 degrees of freedom
## Multiple R-squared:  0.01386,    Adjusted R-squared:  0.01055 
## F-statistic: 4.187 on 1 and 298 DF,  p-value: 0.04162
# Hip
hip_vert_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_hip_vert_LMM)
summary(hip_vert_BA_plot_LR)
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_hip_vert_LMM)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -255.386  -53.520   -0.846   49.427  295.833 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -42.10744   19.93510  -2.112   0.0355 *
## mean          0.03683    0.01685   2.186   0.0296 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 87.33 on 297 degrees of freedom
## Multiple R-squared:  0.01583,    Adjusted R-squared:  0.01252 
## F-statistic: 4.778 on 1 and 297 DF,  p-value: 0.02961

4.3.2 Indices of accuracy

To evaluate models prediction accuracy, MAE, MAPE and RMSE were computed. These indices were computed with the accuracy_indices() function.

# For resultant peak ground reaction force
# Back
back_res_accuracy <- accuracy_indices(LOOCV_back_res_LMM, "pRGRF_N", "pRGRF_N_predicted")
# Hip
hip_res_accuracy <- accuracy_indices(LOOCV_hip_res_LMM, "pRGRF_N", "pRGRF_N_predicted")

# For vertical peak ground reaction force
# Back
back_vert_accuracy <- accuracy_indices(LOOCV_back_vert_LMM, "pVGRF_N", "pVGRF_N_predicted")
# Hip
hip_vert_accuracy <- accuracy_indices(LOOCV_hip_vert_LMM, "pVGRF_N", "pVGRF_N_predicted")

Table below shows these indices for each accelerometer placement in both pRGRF and pVGRF.

accuracy <- data.frame(
  acc_placement = c(rep("Back", 2), rep("Hip", 2)),
  condition     = rep(c("pRGRF", "pVGRF"), 2),
  MAE           = c(
    round(back_res_accuracy[[5]], 1), round(back_vert_accuracy[[5]], 1), 
    round(hip_res_accuracy[[5]], 1), round(hip_vert_accuracy[[5]], 1)
    ),
  MAPE          = c(
    paste(round(back_res_accuracy[[12]], 1), "%", sep = ""), paste(round(back_vert_accuracy[[12]], 1), "%", sep = ""), 
    paste(round(hip_res_accuracy[[12]], 1), "%", sep = ""), paste(round(hip_vert_accuracy[[12]], 1), "%", sep = "")
    ),
  RMSE          = c(
    round(back_res_accuracy[[7]], 1), round(back_vert_accuracy[[7]], 1), 
    round(hip_res_accuracy[[7]], 1), round(hip_vert_accuracy[[7]], 1)
    )
)
kable(accuracy, col.names = c("Accelerometer Placement", "GRF", "MAE", "MAPE", "RMSE"), align = c("l", "l", "r", "r", "r"))
Accelerometer Placement GRF MAE MAPE RMSE
Back pRGRF 66.7 5.9% 88.0
Back pVGRF 74.9 6.5% 101.7
Hip pRGRF 65.5 5.9% 83.6
Hip pVGRF 65.6 5.9% 87.7

All conditions showed a MAE from near 65 N to near 81 N, MAPE from near 6% to near 7% and RMSE from near 84 N to 84 N.

4.3.3 ANOVA

A series of repeated measures analysis of variance (ANOVA) were run to assess whether pGRF predicted from the regression equations were significantly different from those measured with FP. Walking speeds, accelerometer placements (lower back and hip), and the interaction effect (speed X accelerometer placements) were considered for analysis. These procedures were taken separately for resultant and its vertical component. If assumptions of sphericity were violated (p < 0.05), the conservative Greenhouse–Geisser correction factor would be applied to adjust the degrees of freedom. Post-hoc analyses would be conducted using pairwise comparisons with Holm’s correction if a significant difference were observed among actual and predicted pGRF.

In order to run the repeated measures ANOVA, new data frames needed to be built, for both pRGRF and pVGRF, putting all the pGRF values in a single variable and grouping then by accelerometer placement or actual value in another variable. An example of code to tidy the data is shown below.

## Predicted pRGRF
back_res_pred <- LOOCV_back_res_LMM %>% 
  select(ID, speed, pRGRF_N_predicted) %>% 
  spread(key = speed, value = pRGRF_N_predicted) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = back
  )

hip_res_pred <- LOOCV_hip_res_LMM %>% 
  select(ID, speed, pRGRF_N_predicted) %>% 
  spread(key = speed, value = pRGRF_N_predicted) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = hip
  )

res_pred_df <- back_res_pred %>% 
  full_join(hip_res_pred, by = c("ID", "speed")) %>% 
  na.omit()

## Actual pRGRF
back_res_actual <- LOOCV_back_res_LMM %>% 
  select(ID, speed, pRGRF_N) %>% 
  spread(key = speed, value = pRGRF_N) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = actual_back
  )

hip_res_actual <- LOOCV_hip_res_LMM %>% 
  select(ID, speed, pRGRF_N) %>% 
  spread(key = speed, value = pRGRF_N) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = actual_hip
  )

res_actual_df <- back_res_actual %>% 
  full_join(hip_res_actual, by = c("ID", "speed")) %>% 
  na.omit()

res_actual_df <- res_actual_df %>% 
  mutate(actual = (actual_back + actual_hip) / 2) %>% 
  select(ID, speed, actual) 

## Merge predicted and actual data frames
res_ANOVA_df <- res_pred_df %>% 
  full_join(res_actual_df, by = c("ID", "speed")) %>% 
  gather(
    back, hip, actual,
    key = "group",
    value = "pRGRF"
  )
res_ANOVA_df$ID    <- as.factor(res_ANOVA_df$ID)
res_ANOVA_df$speed <- as.factor(res_ANOVA_df$speed)
res_ANOVA_df$group <- as.factor(res_ANOVA_df$group)

Repeated measures ANOVAs were performed using the ezANOVA() function of the ez package. Summary statistics for pRGRF and pVGRF ANOVA can be seen below, respectively.

# For resultant peak ground reaction force
res_ANOVA <- ezANOVA(
  data     = res_ANOVA_df,
  dv       = pRGRF,
  wid      = ID,
  within   = .(speed, group),
  detailed = TRUE,
  type     = 3
)
res_ANOVA
## $ANOVA
##        Effect DFn DFd          SSn        SSd            F             p
## 1 (Intercept)   1  45 8.401226e+08 41784505.3 9.047736e+02  1.922348e-31
## 2       speed   4 180 7.469909e+06   465997.2 7.213474e+02 1.330059e-109
## 3       group   2  90 8.151181e+01   876689.8 4.183956e-03  9.958250e-01
## 4 speed:group   8 360 1.474347e+04   331353.9 2.002259e+00  4.530694e-02
##   p<.05          ges
## 1     * 9.508154e-01
## 2     * 1.466746e-01
## 3       1.875619e-06
## 4     * 3.391386e-04
## 
## $`Mauchly's Test for Sphericity`
##        Effect            W            p p<.05
## 2       speed 0.0353185983 1.104234e-26     *
## 3       group 0.4299955589 8.633907e-09     *
## 4 speed:group 0.0005121196 3.532794e-47     *
## 
## $`Sphericity Corrections`
##        Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF]
## 2       speed 0.4083488 2.563967e-46         * 0.4216157 9.726767e-48
## 3       group 0.6369409 9.746061e-01           0.6471940 9.758931e-01
## 4 speed:group 0.2655135 1.378260e-01           0.2790290 1.350822e-01
##   p[HF]<.05
## 2         *
## 3          
## 4
# For vertical peak ground reaction force
vert_ANOVA <- ezANOVA(
  data     = vert_ANOVA_df,
  dv       = pVGRF,
  wid      = ID,
  within   = .(speed, group),
  detailed = TRUE,
  type     = 3
)
vert_ANOVA
## $ANOVA
##        Effect DFn DFd          SSn        SSd            F             p
## 1 (Intercept)   1  45 8.079216e+08 38787197.4 937.33175604  8.997743e-32
## 2       speed   4 180 5.570015e+06   377459.7 664.04624137 1.443126e-106
## 3       group   2  90 1.426599e+03  1353453.5   0.04743195  9.536992e-01
## 4 speed:group   8 360 7.412484e+03   339268.4   0.98317970  4.487240e-01
##   p<.05          ges
## 1     * 9.518633e-01
## 2     * 1.199726e-01
## 3       3.491533e-05
## 4       1.813905e-04
## 
## $`Mauchly's Test for Sphericity`
##        Effect            W            p p<.05
## 2       speed 0.0449465515 1.555112e-24     *
## 3       group 0.8418464231 2.265287e-02     *
## 4 speed:group 0.0005848457 4.244722e-46     *
## 
## $`Sphericity Corrections`
##        Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF]
## 2       speed 0.4411131 1.738440e-48         * 0.4577530 3.227609e-50
## 3       group 0.8634433 9.346655e-01           0.8947448 9.396637e-01
## 4 speed:group 0.3093066 3.915361e-01           0.3286990 3.951765e-01
##   p[HF]<.05
## 2         *
## 3          
## 4

For either resultant and its vertical component, repeated measures ANOVA demonstrated that actual and predicted pGRF increased significantly (p < 0.001) along with increments in walking speed.


To assess for any significant differences between actual and predicted pGRF in each speed, five separate repeated measures ANOVA were done, one for each walking speed.

# For resultant peak ground reaction force
# 2 km/h
res_ANOVA_s2 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 2),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 3 km/h
res_ANOVA_s3 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 3),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 4 km/h
res_ANOVA_s4 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 4),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 5 km/h
res_ANOVA_s5 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 5),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 6 km/h
res_ANOVA_s6 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 6),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# For vertical peak ground reaction force
# 2 km/h
vert_ANOVA_s2 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 2),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 3 km/h
vert_ANOVA_s3 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 3),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 4 km/h
vert_ANOVA_s4 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 4),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 5 km/h
vert_ANOVA_s5 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 5),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 6 km/h
vert_ANOVA_s6 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 6),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

No significant differences (p > 0.05) between actual and predicted pGRF in each speed were found.

These results can be seen in figure below.

source(here("figs", "fig3.R"))

GRF_plot_grid

4.4 Regressions comparison

Our equation was compared with a previously published reference equation (11), in which Neugebauer et al. used a similar approach for the prediction of pGRF. Comparison was performed using the equation to predict pVGRF from hip-worn accelerometers, as it was the only suitable with our data. The analysis was performed in three ways using: i) a subsample of participants with normal weight or overweight (BMI ≥ 18.5 and < 30 kg.m-2); ii) a subsample of obese participants (BMI ≥ 30 kg.m-2); iii) the whole sample.

New data frames needed to be built for each of the subsamples and the whole sample. Then, pVGRF was predicted using the reference equation. Code to build the data frames is shown below.

# Prepare  data frames
non_obese_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>% 
  select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted) %>% 
  filter(BMI < 30)

obese_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>% 
  select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted) %>% 
  filter(BMI >= 30)

whole_sample_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>% 
  select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted)

# Apply Neugebauer 2014 equation
non_obese_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(non_obese_df)) {
  non_obese_df$pVGRF_N_Neugebauer[i] <-  
    exp(5.247 + (0.271 * non_obese_df$pVACC_g[i]) + (0.014 * non_obese_df$body_mass[i]))
}

obese_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(obese_df)) {
  obese_df$pVGRF_N_Neugebauer[i] <-  
    exp(5.247 + (0.271 * obese_df$pVACC_g[i]) + (0.014 * obese_df$body_mass[i]))
}

whole_sample_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(whole_sample_df)) {
  whole_sample_df$pVGRF_N_Neugebauer[i] <-  
    exp(5.247 + (0.271 * whole_sample_df$pVACC_g[i]) + (0.014 * whole_sample_df$body_mass[i]))
}

Bland-Altman plots were used to confront the agreement between pVGRF measured with FP and those predicted using both regression equations.

source(here("figs", "fig4.R"))

BA_plot_grid

Code to generate these plots can be seen here.

A lower dispersion around the bias value, as well as a lower percentage of values falling off the LoA can be appreciated in our equation. Additionally, while the bias in our equation tended to be zero, in the reference equation bias was always > 0 (p < 0.05), showing a consistent underestimation of pVGRF.


Also, to assess the prediction accuracy, MAE, MAPE and RMSE were calculated using pVGRF values predicted from both equations.

non_obese_our_accuracy     <- accuracy_indices(non_obese_df, "pVGRF_N", "pVGRF_N_predicted")
non_obese_Neug_accuracy    <- accuracy_indices(non_obese_df, "pVGRF_N", "pVGRF_N_Neugebauer")
obese_our_accuracy         <- accuracy_indices(obese_df, "pVGRF_N", "pVGRF_N_predicted")
obese_Neug_accuracy        <- accuracy_indices(obese_df, "pVGRF_N", "pVGRF_N_Neugebauer")
whole_sample_our_accuracy  <- accuracy_indices(whole_sample_df, "pVGRF_N", "pVGRF_N_predicted")
whole_sample_Neug_accuracy <- accuracy_indices(whole_sample_df, "pVGRF_N", "pVGRF_N_Neugebauer")

The accuracy indices from our equation were all substantially lower compared to the reference equation indices, with an overall MAPE approximately 3 times smaller for our equation. For both equations, the MAPE was lower in the obese subsample than in the non-obese subsample

5 Loading rate prediction

5.1 Linear mixed models

For the LR prediction, the same steps were followed as the GRF prediction. The code to generate the models can be seen bellow.

# For resultant peak loading rate
back_res_LR_LMM <- lme(
  fixed = pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass,
  random = ~ 1 | ID,
  method = "ML",
  correlation = corAR1(),
  data = back_LR,
  na.action = na.omit
)
r2_back_res_LR_LMM <- rsquared(back_res_LR_LMM)

summary(back_res_LR_LMM)
## Linear mixed-effects model fit by maximum likelihood
##  Data: back_LR 
##        AIC      BIC   logLik
##   4874.519 4903.336 -2429.26
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.4054883 2372.028
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.6605835 
## Fixed effects: pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass 
##                        Value Std.Error  DF   t-value p-value
## (Intercept)        -287.0209 1590.2580 211 -0.180487  0.8569
## pRATR_gs            572.7967  102.1134 211  5.609416  0.0000
## I(pRATR_gs^2)        -9.8958    1.6666 211 -5.937732  0.0000
## body_mass            18.1178   18.8758  55  0.959841  0.3413
## pRATR_gs:body_mass    3.4078    1.1110 211  3.067340  0.0024
##  Correlation: 
##                    (Intr) pRATR_g I(RATR bdy_ms
## pRATR_gs           -0.720                      
## I(pRATR_gs^2)       0.075 -0.336               
## body_mass          -0.950  0.603   0.146       
## pRATR_gs:body_mass  0.693 -0.810  -0.243 -0.738
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.4916791 -0.5669550 -0.1458794  0.4684454  4.2366891 
## 
## Number of Observations: 271
## Number of Groups: 57
rsquared(back_res_LR_LMM)
##   Response   family     link method  Marginal Conditional
## 1  pRLR_Ns gaussian identity   none 0.7500468   0.7500468
hip_res_LR_LMM <- lme(
  fixed = pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass,
  random = ~ 1 | ID,
  method = "ML",
  correlation = corAR1(),
  data = hip_LR,
  na.action = na.omit
)
r2_hip_res_LR_LMM <- rsquared(hip_res_LR_LMM)

summary(hip_res_LR_LMM)
## Linear mixed-effects model fit by maximum likelihood
##  Data: hip_LR 
##        AIC      BIC    logLik
##   4917.796 4946.584 -2450.898
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.5023566 2529.466
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.6017064 
## Fixed effects: pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass 
##                        Value Std.Error  DF   t-value p-value
## (Intercept)        -3510.410 1846.9416 210 -1.900661  0.0587
## pRATR_gs             514.898  110.2864 210  4.668734  0.0000
## I(pRATR_gs^2)         -8.639    1.4312 210 -6.035804  0.0000
## body_mass             51.937   20.6749  55  2.512072  0.0150
## pRATR_gs:body_mass     2.929    1.0337 210  2.833284  0.0051
##  Correlation: 
##                    (Intr) pRATR_g I(RATR bdy_ms
## pRATR_gs           -0.785                      
## I(pRATR_gs^2)       0.329 -0.615               
## body_mass          -0.948  0.657  -0.096       
## pRATR_gs:body_mass  0.748 -0.829   0.104 -0.780
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.22341333 -0.61951158 -0.07192759  0.43894165  3.05373725 
## 
## Number of Observations: 270
## Number of Groups: 57
rsquared(hip_res_LR_LMM)
##   Response   family     link method  Marginal Conditional
## 1  pRLR_Ns gaussian identity   none 0.7021441   0.7021441
# For vertical peak loading rate
back_vert_LR_LMM <- lme(
  fixed = pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass,
  random = ~ 1 | ID,
  method = "ML",
  correlation = corAR1(),
  data = back_LR,
  na.action = na.omit
)
r2_back_vert_LR_LMM <- rsquared(back_vert_LR_LMM)

summary(back_vert_LR_LMM)
## Linear mixed-effects model fit by maximum likelihood
##  Data: back_LR 
##        AIC      BIC    logLik
##   4858.839 4887.626 -2421.419
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.4732634 2393.903
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.6661911 
## Fixed effects: pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass 
##                        Value Std.Error  DF   t-value p-value
## (Intercept)        -324.0761 1622.7623 210 -0.199706  0.8419
## pVATR_gs            552.8242  103.4180 210  5.345534  0.0000
## I(pVATR_gs^2)       -11.9453    1.8892 210 -6.322881  0.0000
## body_mass            18.1405   19.3620  55  0.936915  0.3529
## pVATR_gs:body_mass    3.9586    1.1637 210  3.401656  0.0008
##  Correlation: 
##                    (Intr) pVATR_g I(VATR bdy_ms
## pVATR_gs           -0.722                      
## I(pVATR_gs^2)       0.031 -0.287               
## body_mass          -0.949  0.604   0.193       
## pVATR_gs:body_mass  0.689 -0.790  -0.326 -0.743
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.07632135 -0.62184789 -0.08981213  0.46058800  3.09934690 
## 
## Number of Observations: 270
## Number of Groups: 57
rsquared(back_vert_LR_LMM)
##   Response   family     link method  Marginal Conditional
## 1  pVLR_Ns gaussian identity   none 0.7283041   0.7283041
hip_vert_LR_LMM <- lme(
  fixed = pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass,
  random = ~ 1 | ID,
  method = "ML",
  correlation = corAR1(),
  data = hip_LR,
  na.action = na.omit
)
r2_hip_vert_LR_LMM <- rsquared(hip_vert_LR_LMM)

summary(hip_vert_LR_LMM)
## Linear mixed-effects model fit by maximum likelihood
##  Data: hip_LR 
##        AIC      BIC    logLik
##   4907.858 4936.645 -2445.929
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.4267048 2637.244
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.6724967 
## Fixed effects: pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass 
##                         Value Std.Error  DF   t-value p-value
## (Intercept)        -2687.8662 1930.3344 210 -1.392436  0.1653
## pVATR_gs             407.8434  112.4430 210  3.627111  0.0004
## I(pVATR_gs^2)         -7.6603    1.4186 210 -5.400053  0.0000
## body_mass             45.8905   21.7542  55  2.109503  0.0395
## pVATR_gs:body_mass     3.8995    1.0739 210  3.631057  0.0004
##  Correlation: 
##                    (Intr) pVATR_g I(VATR bdy_ms
## pVATR_gs           -0.763                      
## I(pVATR_gs^2)       0.305 -0.596               
## body_mass          -0.951  0.643  -0.086       
## pVATR_gs:body_mass  0.729 -0.839   0.101 -0.756
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.4906751 -0.5800428 -0.1047851  0.4384381  2.9373647 
## 
## Number of Observations: 270
## Number of Groups: 57
rsquared(hip_vert_LR_LMM)
##   Response   family     link method  Marginal Conditional
## 1  pVLR_Ns gaussian identity   none 0.6868396   0.6868396

In all models, conditional R2 values ranged from 0.69 to 0.75.

5.2 Validation analysis

Similar to the GRF prediction models, LR model validation was assessed by the LOOCV method with the cross_validate_mixed_model() function and is shown below.

# For resultant peak loading rate
fix_eff    <- pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass
rand_eff   <- ~ 1 | ID
# Back
LOOCV_back_res_LR_LMM <- do.call(rbind, (lapply(unique(back_LR$ID), cross_validate_mixed_model, df = back_LR)))
# Hip
LOOCV_hip_res_LR_LMM <- do.call(rbind, (lapply(unique(hip_LR$ID), cross_validate_mixed_model, df = hip_LR)))

# For vertical peak loading rate
fix_eff    <- pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass
rand_eff   <- ~ 1 | ID
# Back
LOOCV_back_vert_LR_LMM <- do.call(rbind, (lapply(unique(back_LR$ID), cross_validate_mixed_model, df = back_LR)))
# Hip
LOOCV_hip_vert_LR_LMM <- do.call(rbind, (lapply(unique(hip_LR$ID), cross_validate_mixed_model, df = hip_LR)))

5.2.1 Bland-Altman plots

Bland-Altman plots were also used to examine the agreement between actual and predicted pLR for both resultant and its vertical component.

Bland-Altman plots can be seen below. Panel C shows plots of pRLR and panel D shows plots of pVLR.

source(here("figs", "fig2.R"))

BA_plot_grid_1 <- plot_grid(
  back_res_LR_BA_plot  + theme(legend.position = "none"),
  back_vert_LR_BA_plot + theme(legend.position = "none"),
  hip_res_LR_BA_plot   + theme(legend.position = "none"),
  hip_vert_LR_BA_plot  + theme(legend.position = "none"),
  labels = c("C", "D", "", ""),
  align  = "h", vjust = 1, label_size = 16,
  ncol   = 2, nrow = 2
)

legend <- get_legend(back_res_LR_BA_plot)

BA_plot_grid <- plot_grid(BA_plot_grid_1, legend, ncol = 1, rel_heights = c(1, 0.1))

BA_plot_grid

Code to generate these plots can be seen here.


One-sample t-tests were run to check whether bias from each accelerometer placement, in both pRGRF and pVGRF, was significantly different than 0. These tests were run using the t.test() function of the base R with the argument mu = 0. No significant differences were found.

## 
##  One Sample t-test
## 
## data:  LOOCV_back_res_LR_LMM$diff
## t = -0.95034, df = 269, p-value = 0.3428
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -418.0833  145.8665
## sample estimates:
## mean of x 
## -136.1084
## 
##  One Sample t-test
## 
## data:  LOOCV_hip_res_LR_LMM$diff
## t = -0.73462, df = 269, p-value = 0.4632
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -418.3023  190.9683
## sample estimates:
## mean of x 
##  -113.667
## 
##  One Sample t-test
## 
## data:  LOOCV_back_vert_LR_LMM$diff
## t = -0.97483, df = 269, p-value = 0.3305
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -428.7289  144.7705
## sample estimates:
## mean of x 
## -141.9792
## 
##  One Sample t-test
## 
## data:  LOOCV_hip_vert_LR_LMM$diff
## t = -1.1318, df = 269, p-value = 0.2587
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -497.1642  134.2107
## sample estimates:
## mean of x 
## -181.4768

Also, linear regressions were applied to identify if there was any proportional bias, that is, if bias was related with the magnitude of the mean between measured and predicted pGRF (10). Linear regressions were run using the lm() function of the base R.

Results showed a significant proportional bias (p < 0.05), however, with a low magnitude. These results suggest that despite there is a trend for underestimation at increasingly higher pLR values, the magnitude of this effect is neglectable.

# For resultant peak loading rate
# Back
back_res_LR_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_res_LR_LMM)
# Hip
hip_res_LR_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_hip_res_LR_LMM)

# For vertical peak loading rate
# Back
back_vert_LR_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_vert_LR_LMM)
# Hip
hip_vert_LR_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_hip_vert_LR_LMM)

5.2.2 Indices of accuracy

To evaluate models prediction accuracy, MAE, MAPE and RMSE were computed. These indices were computed with the accuracy_indices() function.

# For resultant peak loading rate
# Back
back_res_LR_accuracy <- accuracy_indices(LOOCV_back_res_LR_LMM, "pRLR_Ns", "pRLR_Ns_predicted")
# Hip
hip_res_LR_accuracy <- accuracy_indices(LOOCV_hip_res_LR_LMM, "pRLR_Ns", "pRLR_Ns_predicted")

# For vertical peak loading rate
# Back
back_vert_LR_accuracy <- accuracy_indices(LOOCV_back_vert_LR_LMM, "pVLR_Ns", "pVLR_Ns_predicted")
# Back
hip_vert_LR_accuracy <- accuracy_indices(LOOCV_hip_vert_LR_LMM, "pVLR_Ns", "pVLR_Ns_predicted")

Table below shows these indices for each accelerometer placement in both pRLR and pVLR.

accuracy <- data.frame(
  acc_placement = c(rep("Back", 2), rep("Hip", 2)),
  condition     = rep(c("pRLR", "pVLR"), 2),
  MAE           = c(
    round(back_res_LR_accuracy[[5]], 1), round(back_vert_LR_accuracy[[5]], 1), 
    round(hip_res_LR_accuracy[[5]], 1), round(hip_vert_LR_accuracy[[5]], 1)
    ),
  MAPE          = c(
    paste(round(back_res_LR_accuracy[[12]], 1), "%", sep = ""), paste(round(back_vert_LR_accuracy[[12]], 1), "%", sep = ""), 
    paste(round(hip_res_LR_accuracy[[12]], 1), "%", sep = ""), paste(round(hip_vert_LR_accuracy[[12]], 1), "%", sep = "")
    ),
  RMSE          = c(
    round(back_res_LR_accuracy[[7]], 1), round(back_vert_LR_accuracy[[7]], 1), 
    round(hip_res_LR_accuracy[[7]], 1), round(hip_vert_LR_accuracy[[7]], 1)
    )
)
kable(accuracy, col.names = c("Accelerometer Placement", "LR", "MAE", "MAPE", "RMSE"), align = c("l", "l", "r", "r", "r"))
Accelerometer Placement LR MAE MAPE RMSE
Back pRLR 1732.7 19.5% 2352.9
Back pVLR 1771.9 19.8% 2393.0
Hip pRLR 1866.8 20.6% 2540.3
Hip pVLR 1975.5 22.4% 2636.1

5.2.3 ANOVA

Also, to assess if pLR predicted values significantly differed from the actual values measured by the FP, repeated measures ANOVA were run. As with the ANOVA for the pGRF, here walking speeds, accelerometer placements (lower back and hip), and the interaction effect (speed X accelerometer placements) were considered for analysis.

Code to tidy the data is shown below.

## Predicted pRLR
back_res_LR_pred <- LOOCV_back_res_LR %>% 
  select(ID, speed, pRLR_Ns_predicted) %>% 
  spread(key = speed, value = pRLR_Ns_predicted) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = back
  )

hip_res_LR_pred <- LOOCV_hip_res_LR %>%
  select(ID, speed, pRLR_Ns_predicted) %>% 
  spread(key = speed, value = pRLR_Ns_predicted) %>%
  na.omit() %>%
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = hip
  )

res_pred_LR_df <- back_res_LR_pred %>% 
  full_join(hip_res_LR_pred, by = c("ID", "speed")) %>% 
  na.omit()

## Actual pRLR
back_res_actual_LR <- LOOCV_back_res_LR %>% 
  select(ID, speed, pRLR_Ns) %>% 
  spread(key = speed, value = pRLR_Ns) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = actual_back
  )

hip_res_actual_LR <- LOOCV_hip_res_LR %>% 
  select(ID, speed, pRLR_Ns) %>% 
  spread(key = speed, value = pRLR_Ns) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = actual_hip
  )

res_actual_LR_df <- back_res_actual_LR %>% 
  full_join(hip_res_actual_LR, by = c("ID", "speed")) %>% 
  na.omit()

res_actual_LR_df <- res_actual_LR_df %>% 
  mutate(actual = (actual_back + actual_hip) / 2) %>% 
  select(ID, speed, actual) 

## Merge predicted and actual data frames
res_ANOVA_LR_df <- res_pred_LR_df %>% 
  full_join(res_actual_LR_df, by = c("ID", "speed")) %>% 
  gather(
    back, hip, actual,
    key = "group",
    value = "pRLR"
  )
res_ANOVA_LR_df$ID    <- as.factor(res_ANOVA_LR_df$ID)
res_ANOVA_LR_df$speed <- as.factor(res_ANOVA_LR_df$speed)
res_ANOVA_LR_df$group <- as.factor(res_ANOVA_LR_df$group)

Repeated measures ANOVAs were performed using the ezANOVA() function of the ez package. Code to run the analysis is shown below.

## 
##  Pairwise comparisons using paired t tests 
## 
## data:  s2_res$pRLR and s2_res$group 
## 
##      actual  back   
## back 0.00016 -      
## hip  0.00011 0.72304
## 
## P value adjustment method: holm

These results can be seen in figure below.

source(here("figs", "fig3.R"))

LR_plot_grid

6 R session info

devtools::session_info()
## ─ Session info ──────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.5.2 (2018-12-20)
##  os       macOS  10.15.1              
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Europe/Lisbon               
##  date     2019-12-23                  
## 
## ─ Packages ──────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  abind          1.4-5    2016-07-21 [1] CRAN (R 3.5.0)
##  acepack        1.4.1    2016-10-29 [1] CRAN (R 3.5.0)
##  assertthat     0.2.0    2017-04-11 [1] CRAN (R 3.5.0)
##  backports      1.1.2    2017-12-13 [1] CRAN (R 3.5.0)
##  base64enc      0.1-3    2015-07-28 [1] CRAN (R 3.5.0)
##  broman       * 0.68-2   2018-07-25 [1] CRAN (R 3.5.0)
##  broom          0.4.5    2018-07-03 [1] CRAN (R 3.5.0)
##  callr          2.0.4    2018-05-15 [1] CRAN (R 3.5.0)
##  car            3.0-2    2018-08-23 [1] CRAN (R 3.5.0)
##  carData        3.0-2    2018-09-30 [1] CRAN (R 3.5.0)
##  cellranger     1.1.0    2016-07-27 [1] CRAN (R 3.5.0)
##  checkmate      1.8.5    2017-10-24 [1] CRAN (R 3.5.0)
##  cli            1.0.1    2018-09-25 [1] CRAN (R 3.5.0)
##  cluster        2.0.7-1  2018-04-13 [1] CRAN (R 3.5.2)
##  colorspace     1.3-2    2016-12-14 [1] CRAN (R 3.5.0)
##  cowplot      * 0.9.3    2018-07-15 [1] CRAN (R 3.5.0)
##  crayon         1.3.4    2017-09-16 [1] CRAN (R 3.5.0)
##  curl           3.2      2018-03-28 [1] CRAN (R 3.5.0)
##  data.table     1.12.0   2019-01-13 [1] CRAN (R 3.5.2)
##  desc           1.2.0    2018-05-01 [1] CRAN (R 3.5.0)
##  devtools       2.0.1    2018-10-26 [1] CRAN (R 3.5.0)
##  digest         0.6.19   2019-05-20 [1] CRAN (R 3.5.2)
##  dplyr        * 0.8.0.1  2019-02-15 [1] CRAN (R 3.5.2)
##  evaluate       0.10.1   2017-06-24 [1] CRAN (R 3.5.0)
##  ez           * 4.4-0    2016-11-02 [1] CRAN (R 3.5.0)
##  fansi          0.4.0    2018-10-05 [1] CRAN (R 3.5.0)
##  forcats      * 0.3.0    2018-02-19 [1] CRAN (R 3.5.0)
##  foreign        0.8-71   2018-07-20 [1] CRAN (R 3.5.2)
##  Formula        1.2-3    2018-05-03 [1] CRAN (R 3.5.0)
##  fs             1.2.6    2018-08-23 [1] CRAN (R 3.5.0)
##  ggplot2      * 3.1.0    2018-10-25 [1] CRAN (R 3.5.0)
##  glue           1.3.0    2018-07-17 [1] CRAN (R 3.5.0)
##  gridExtra      2.3      2017-09-09 [1] CRAN (R 3.5.0)
##  gtable         0.2.0    2016-02-26 [1] CRAN (R 3.5.0)
##  haven          2.0.0    2018-11-22 [1] CRAN (R 3.5.0)
##  here         * 0.1      2017-05-28 [1] CRAN (R 3.5.0)
##  highr          0.7      2018-06-09 [1] CRAN (R 3.5.0)
##  Hmisc          4.3-0    2019-11-07 [1] CRAN (R 3.5.2)
##  hms            0.4.2    2018-03-10 [1] CRAN (R 3.5.0)
##  htmlTable      1.12     2018-05-26 [1] CRAN (R 3.5.0)
##  htmltools      0.3.6    2017-04-28 [1] CRAN (R 3.5.0)
##  htmlwidgets    1.3      2018-09-30 [1] CRAN (R 3.5.0)
##  httr           1.3.1    2017-08-20 [1] CRAN (R 3.5.0)
##  jsonlite       1.5      2017-06-01 [1] CRAN (R 3.5.0)
##  knitr        * 1.21     2018-12-10 [1] CRAN (R 3.5.0)
##  labeling       0.3      2014-08-23 [1] CRAN (R 3.5.0)
##  lattice        0.20-38  2018-11-04 [1] CRAN (R 3.5.2)
##  latticeExtra   0.6-28   2016-02-09 [1] CRAN (R 3.5.0)
##  lazyeval       0.2.1    2017-10-29 [1] CRAN (R 3.5.0)
##  lme4           1.1-19   2018-11-10 [1] CRAN (R 3.5.0)
##  lubridate      1.7.4    2018-04-11 [1] CRAN (R 3.5.0)
##  magrittr       1.5      2014-11-22 [1] CRAN (R 3.5.0)
##  MASS           7.3-51.1 2018-11-01 [1] CRAN (R 3.5.2)
##  Matrix         1.2-15   2018-11-01 [1] CRAN (R 3.5.2)
##  memoise        1.1.0    2017-04-21 [1] CRAN (R 3.5.0)
##  mgcv           1.8-26   2018-11-21 [1] CRAN (R 3.5.2)
##  minqa          1.2.4    2014-10-09 [1] CRAN (R 3.5.0)
##  mnormt         1.5-5    2016-10-15 [1] CRAN (R 3.5.0)
##  modelr         0.1.2    2018-05-11 [1] CRAN (R 3.5.0)
##  munsell        0.5.0    2018-06-12 [1] CRAN (R 3.5.0)
##  nlme         * 3.1-137  2018-04-07 [1] CRAN (R 3.5.2)
##  nloptr         1.2.1    2018-10-03 [1] CRAN (R 3.5.0)
##  nnet           7.3-12   2016-02-02 [1] CRAN (R 3.5.0)
##  openxlsx       4.1.0.1  2019-05-28 [1] CRAN (R 3.5.2)
##  pbkrtest       0.4-7    2017-03-15 [1] CRAN (R 3.5.0)
##  piecewiseSEM * 2.0.2    2018-07-24 [1] CRAN (R 3.5.0)
##  pillar         1.3.1    2018-12-15 [1] CRAN (R 3.5.0)
##  pkgbuild       1.0.2    2018-10-16 [1] CRAN (R 3.5.0)
##  pkgconfig      2.0.2    2018-08-16 [1] CRAN (R 3.5.0)
##  pkgload        1.0.2    2018-10-29 [1] CRAN (R 3.5.0)
##  plyr           1.8.4    2016-06-08 [1] CRAN (R 3.5.0)
##  prettyunits    1.0.2    2015-07-13 [1] CRAN (R 3.5.0)
##  processx       3.1.0    2018-05-15 [1] CRAN (R 3.5.0)
##  psych          1.8.4    2018-05-06 [1] CRAN (R 3.5.0)
##  purrr        * 0.3.0    2019-01-27 [1] CRAN (R 3.5.2)
##  R6             2.4.0    2019-02-14 [1] CRAN (R 3.5.2)
##  RColorBrewer   1.1-2    2014-12-07 [1] CRAN (R 3.5.0)
##  Rcpp           1.0.0    2018-11-07 [1] CRAN (R 3.5.0)
##  readr        * 1.3.1    2018-12-21 [1] CRAN (R 3.5.0)
##  readxl         1.1.0    2018-04-20 [1] CRAN (R 3.5.0)
##  remotes        2.0.2    2018-10-30 [1] CRAN (R 3.5.0)
##  reshape2       1.4.3    2017-12-11 [1] CRAN (R 3.5.0)
##  rio            0.5.16   2018-11-26 [1] CRAN (R 3.5.0)
##  rlang          0.4.1    2019-10-24 [1] CRAN (R 3.5.2)
##  rmarkdown      1.10     2018-06-11 [1] CRAN (R 3.5.0)
##  rpart          4.1-13   2018-02-23 [1] CRAN (R 3.5.2)
##  rprojroot      1.3-2    2018-01-03 [1] CRAN (R 3.5.0)
##  rstudioapi     0.7      2017-09-07 [1] CRAN (R 3.5.0)
##  rvest          0.3.2    2016-06-17 [1] CRAN (R 3.5.0)
##  scales         1.0.0    2018-08-09 [1] CRAN (R 3.5.0)
##  sessioninfo    1.1.1    2018-11-05 [1] CRAN (R 3.5.0)
##  stringi        1.3.1    2019-02-13 [1] CRAN (R 3.5.2)
##  stringr      * 1.4.0    2019-02-10 [1] CRAN (R 3.5.2)
##  survival       2.43-3   2018-11-26 [1] CRAN (R 3.5.2)
##  testthat       2.0.0    2017-12-13 [1] CRAN (R 3.5.0)
##  tibble       * 2.0.1    2019-01-12 [1] CRAN (R 3.5.2)
##  tidyr        * 0.8.2    2018-10-28 [1] CRAN (R 3.5.0)
##  tidyselect     0.2.5    2018-10-11 [1] CRAN (R 3.5.0)
##  tidyverse    * 1.2.1    2017-11-14 [1] CRAN (R 3.5.0)
##  usethis        1.4.0    2018-08-14 [1] CRAN (R 3.5.0)
##  utf8           1.1.4    2018-05-24 [1] CRAN (R 3.5.0)
##  withr          2.1.2    2018-03-15 [1] CRAN (R 3.5.0)
##  xfun           0.4      2018-10-23 [1] CRAN (R 3.5.0)
##  xml2           1.2.0    2018-01-24 [1] CRAN (R 3.5.0)
##  yaml           2.1.19   2018-05-01 [1] CRAN (R 3.5.0)
##  zip            1.0.0    2017-04-25 [1] CRAN (R 3.5.0)
## 
## [1] /Library/Frameworks/R.framework/Versions/3.5/Resources/library

References

1. Guthold R, Stevens GA, Riley LM, Bull FC. Worldwide trends in insufficient physical activity from 2001 to 2016: a pooled analysis of 358 population-based surveys with 19 million participants. Lancet Glob Heal. 2018;6(10):e1077–86.

2. Li X, Gong X, Jiang W. Abdominal obesity and risk of hip fracture: a meta-analysis of prospective studies. Osteoporos Int. 2017;28(10):2747–57.

3. Hind K, Rudman H, Pearce M, Treadgold L, Birrell F. Obesity, bone density for weight and prevalent vertebral fracture at age 62 years: the Newcastle Thousand Families Study. J Clin Densitom. 2018;21(4):606.

4. Kaze AD, Rosen HN, Paik JM. A meta-analysis of the association between body mass index and risk of vertebral fracture. Osteoporos Int. 2018;29(1):31–9.

5. Viester L, Verhagen EA, Hengel KM, Koppes LL, Van Der Beek AJ, Bongers PM. The relation between body mass index and musculoskeletal symptoms in the working population. BMC Musculoskelet Disord. 2013;14(238).

6. Nakagawa S, Schielzeth H, O’Hara RB. A general and simple method for obtainingR2from generalized linear mixed-effects models. Methods in Ecology and Evolution. 2013;4(2):133–42.

7. Jaeger BC, Edwards LJ, Das K, Sen PK. An r2 statistic for fixed effects in the generalized linear mixed model. Journal of Applied Statistics. 2016;44(6):1086–105.

8. Staudenmayer J, Zhu W, Catellier DJ. Statistical considerations in the analysis of accelerometry-based activity monitor data. Medicine & Science in Sports & Exercise. 2012;44(1 Suppl 1):S61–7.

9. Bland JM, Altman DG. Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet. 1986;1(8476):301–10.

10. Giavarina D. Understanding bland altman analysis. Biochemia Medica. 2015;25(2):141–51.

11. Neugebauer JM, Collins KH, Hawkins DA. Ground reaction force estimates from ActiGraph GT3X+ hip accelerations. PLoS One. 2014/06/11. 2014;9(6):e99023.